Broyden's Method in Nuclear Structure Calculations 
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Broyden's method, widely used in quantum chemistry electronic-structure calculations for the 
numerical solution of nonlinear equations in many variables, is applied in the context of the nuclear 
many-body problem. Examples include the unitary gas problem, the nuclear density functional 
theory with Skyrme functionals, and the nuclear coupled-cluster theory. The stability of the method, 
its ease of use, and its rapid convergence rates make Broyden's method a tool of choice for large-scale 
nuclear structure calculations. 
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I. INTRODUCTION 

The nuclear many-body problem is undergoing a re- 
naissance. Hand in hand with the experimental develop- 
ments in the science of rare isotopes, a qualitative change 
in theoretical modeling of the nucleus is taking place. De- 
velopments of powerful conceptual, analytic, algorithmic, 
and computational tools enable scientists to probe the 
inner workings of nuclei with far greater precision than 
previously possible. These new tools, including terascale 
computer platforms, make researchers optimistic that the 
goal of developing a comprehensive, quantitative, and 
predictive theory of the nucleus and nucleonic matter is 
achievable 0, 0|- 

The research challenges faced by the nuclear structure 
community are often interdisciplinary. Indeed, the quan- 
tum many-body problem represents one of the great intel- 
lectual and numerical challenges for nuclear and hadronic 
physics, quantum chemistry, atomic and condensed mat- 
ter physics, and materials sciences. Often, to solve prob- 
lems in many-body science, close interactions of physi- 
cists with computer scientists and mathematicians are 
required. 

A theoretical framework aiming at the microscopic de- 
scription of many-body systems and capable of extrapo- 
lating into unknown regions must fulfill several require- 
ments. Namely, (i) it must be general enough to be con- 
fidently applied to unknown species whose properties are 
largely unknown; (ii) it should be capable of handling 
symmetry-breaking effects inherent to finite systems; (iii) 
it should be able to describe both finite systems and bulk 
matter; and (iv) in addition to observables, the method 
should provide associated error bars. These requirements 
are met by the Density Functional Theory (dft) in the 
formulation of Kohn and Sham Q . 

The generalization of the dft to the case of fermionic 



pairing was formulated for electronic superconductors [J] ■ 
The resulting Hartree-Fock-Bogoliubov (hfb) or Bogoli- 
ubov de-Gennes (Bdc) equations can be viewed as the 
generalized Kohn-Sham equations of the standard dft. 

The nuclear dft provides a reliable method for calcu- 
lating properties of nuclei across the whole nuclear mass 
chart. The main ingredient of the nuclear dft Q is the 
energy density functional that depends on proton and 
neutron densities and currents representing distributions 
of nucleonic matter, spins, momentum, and kinetic en- 
ergy, as well as their derivatives (gradient terms). If pair- 
ing correlations are present, the particle densities are aug- 
mented by pairing densities representing correlated nu- 
cleonic pairs Initially, attempts to build effective en- 
ergy functionals were rooted in the empirical zero-range 
Skyrme interaction treated within the Hartree-Fock (hf) 
or HFB approximation. Following the success of dft in 
atomic and molecular physics and, more recently, in the 
physics of ultra-cold atoms, it was realized that the inter- 
action could be secondary to the functional, and recent 
progress has been directed towards extensions of the en- 
ergy functionals to less traditional forms. 

In parallel, much effort is put into trying to connect 
the nuclear dft to more fundamental many-body theory 
based on inter-nucleon interactions. In order to extend 
the ab-initio program to the medium-mass region of the 
nuclear chart, a method which scales softly with system 
size and has a controllable error estimate is called for. 
The Coupled-Cluster (cc) theory is ideal in this respect. 
It scales polynomial with system size and it is size ex- 
tensive, meaning that the energy scales correctly with 
increasing system size. 

This work represents a collaborative effort performed 
under the Universal Nuclear Energy Density Functional 
(unedf) project The goal of unedf is to develop a 
new-generation theoretical framework that will describe 
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all nuclei and their reactions by making use of the massive 
computer power. One of the objectives of unedf is to 
develop computational infrastructure for leadership-class 
computers; hence the algorithmic and computational fo- 
cus of this paper. 

The common feature of dft and CC is to solve the 
quantum many-body problem for interacting fermions 
through self-consistent convergence schemes. The com- 
putational cost of these iterative procedures can become 
very expensive, especially when the size of the model 
space or the number of nuclei processed simultaneously 
is huge. The advent of teraflop computers make such 
large-scale calculations feasible, but in order to take full 
advantage of unique resources, better convergence algo- 
rithms are called for. The aim of the present paper is to 
report the implementation of Broyden's method y, ITo| , a 
quasi-Newton algorithm to solve large sets of non-linear 
equations, in nuclear structure calculations. The hfb (or 
Bdc) equations - the Kohn-Sham equations of the su- 
perfluid dft - represent a coupled system of non-linear 
equations for densities or mean fields. The non-linearity 
enters through the self-consistency, i.e., the dependence 
of mean fields on densities. The CC equations have in fact 
a very similar structure. Here, the non-linearity is hidden 
in the CC intermediates which depend on the particle- hole 
amplitudes. 

Self-consistent formulations of these approaches usu- 
ally diverge when using straight iterations. The simplest 
method commonly used to avoid divergences is the so- 
called linear mixing, in which input and output at a given 
iteration of the self-consistent process are mixed to pro- 
vide input for the next step. We show that Broyden's 
method leads to much improved convergence rates. In 
many cases, the number of iterations required to reach 
convergence is lower by a factor of 3 to 10 with respect 
to linear-mixing iterative schemes. 

The paper is organized as follows. Broyden's method is 
outlined in Sec.|TTl both in its standard version fSec. lII A]) 
that applies to a small number of unknowns, and in a 
modified version (Sec. IIIB|) which is used to crunch very 
large sets of nonlinear equations. The examples of imple- 
mentations come next. Section IIIII discusses the perfor- 
mance of the standard Broyden's method for the case of 
the unitary gas of two fermionic superfluids in spherical 
geometry. In the following examples the modified Broy- 
den's method has been used, as the size of the problem 
is intractable with the standard approach. The general 
case of the symmetry-unrestricted Skyrme-DFT theory 
is discussed in Sec. IIVI Section |V] illustrates the par- 
ticular case of an axially-symmetric Skyrme-DFT prob- 
lem and discusses the ability of the method to converge 
to the local energy extrema (i.e, minima, maxima, and 
saddle points). The performance of the method in ab- 
initio coupled-cluster calculations, where the number of 
unknowns can be as large as 10^, is presented in Sec. IVII 
Finally, the conclusions of our work are contained in 
Sec. ED 



II. BROYDEN'S MIXING METHOD 

When solving a system of self-consistent equations, one 
usually begins with a set of initial conditions (in nuclear 
structure: single-particle wave functions, densities, fields, 
etc.) that linearize the problem. These initial conditions 
can be represented formally by an A^-dimensional vector 
Vj^\ Solving the self-consistent equations with these 
initial conditions defines a new vector 

yj:;) = i{vt^) (1) 

where I{V) is a function of the initial conditions. The 
self-consistency condition requires that the solution V* 
be a fixed-point of this iteration: I{V*) — V* . 

The goal of a fixed-point iteration algorithm is to con- 
struct a new guess 'V^,^™'''^^ in such a way that the itera- 
tion converges, 

F(")^yi™)--Kf)^o, (2) 

within the required tolerance. Using the previous output 
as the next input, V^^™'*'^^ — v}^\ does not guarantee 
convergence. If the resulting divergence exhibits oscilla- 
tory behavior, the usual ansatz is to mix the input at 
iteration m to define the input at iteration m + l: 

V^^'^ ^ aV^^:^^ + {I - a)Vt^ (3a) 
= V^^^+aF^"'K (3b) 

This simple mixing slows down the iteration, and with a 
suitable choice of the mixing parameter a e [0, 1], con- 
vergence may be achieved. However, even if one can find 
a that leads to a convergent solution, there are many in- 
stances where the convergence is very slow. In large-scale 
calculations involving huge numbers of cases, the running 
time to find poorly converging solutions may become pro- 
hibitive. 

If the vector field I{V) is differentiable, one can use 
derivative information to improve convergence. The idea 
is to regard the self-consistent condition — as a 

set of non- linear equations, and to approach the problem 
from the point of view of finding the root: i^(V*) = 
I{V*) — V* = 0. This set of non- linear equations can 
be solved using the multidimensional Newton-Raphson 
method, in which the next guess can be approximated 
by 

yim+l) ^ y{m} _ 

where B^™) = (j(™))-i and jj"^ = dFj"'\v)/dVk is 
the Jacobian matrix of the non-linear equations at V — 
1^1^"*''. For sufficiently smooth functions, this method has 
a quadratic convergence. However, this approach is very 
expensive as it needs the explicit derivative evaluation. 
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A. Standard Broyden's Method 



with m = max(l, m — M) and 



In many situations, the quick evaluation of the inverse 
Jacobian is not possible. In such cases one can often 
achieve superlinear convergence by using a multidimen- 
sional generalization of the secant method, whereby an 
approximate matrix B^™^ is computed using a secant ap- 
proximation based on the information obtained in previ- 
ous iterations. This method is usually referred to as the 
quasi-Newton method. The difficulty in higher dimen- 
sions is that the secant update is not unique, and naive 
implementations fail to converge. Broyden [3] suggested 
using the Sherman-Morrison formula to form an update 
to the inverse of the Jacobian that has good convergence 
properties: 

where 

\SF) = - F*™^. 

One can start the iterative process with the initial guess 
B(°) = al, as this is equivalent to the simple mixing 
At any given step, one may also check the full quasi- 
Newton step (|4]): if it produces larger residuals 
or increases the energy of the solution, then one might 
revert to simple mixing ([3]). This can be especially im- 
portant if the function F{V) varies significantly with the 
step size because the secant approximation may then be 
extremely poor. More sophisticated, but globally conver- 
gent algorithms, are discussed in 



B. Modified Broyden's Method 

If the self-consistent equations involve many variables, 
i.e., N is large, storing the N x N matrix elements of the 
inverse Jacobian and performing the N xN matrix multi- 
plications might be prohibitive; hence, further improve- 
ments are needed. To deal with the size issues and to 
improve performance, Broyden's method has been mod- 
ified 0, [l^, The modification by Johnson is 
now widely used in quantum chemistry. This modified 
Broyden's method only utilizes information obtained in 
M previous iterations, and this information is used in the 
update of the approximate inverse Jacobian. We recall 
here the final expressions of this modified Broyden mix- 
ing procedure (a detailed derivation is given in Ref. p^): 



where 



m — 1 
k—rh 

= aAF(") -I- AF(") 



(7) 



yin+l) _ yin) 



p{n+l) 



_P'(n) 



(8) 



The weights w„ (n = 1, ...,M) are associated with each 
previous iteration and the values Wn = 1 are usually cho- 
sen. The weight wq is assigned to the error in the inverse 
Jacobian and the value wo = 0.01 proposed in gives 
stable results. The first two terms in Eq. ^ are sim- 
ply the linear mixing of Eq. ^ with a mixing parameter 
a, while the last term is an additional correction. As 
mentioned in the parameter a can be chosen to be 
rather large (in the examples of this paper we have used 
a — 0.7) compared to the usual values used in the simple 
mixing approach. 

Equations are all that is required to update V. 

In the modified Broyden's method, large N x N matrices 
do not appear. It is only one M x M matrix a^; and M 
vectors of length N that need to be stored. The iterative 
procedure starts with an initial input guess VJ^*^^ which 
generates the first output solution v}^2- the next it- 
eration, m=l, VilP is a linear combination of and 
v}^} and the correction term is zero. For m=2, the new 

(2) 

V^^ contains the Broyden correction including the infor- 
mation obtained in the previous (m=l) iterations. The 
process is repeated until convergence is reached. 

Our actual implementation of the modified Broy- 
den's method is based on a modified subroutine from 
the package pwscf (Plane-Wave Self-Consistent Field) 
for electronic-structure dft calculations using pseudo- 
potentials and a plane- wave basis set [l^j. The code is 
published under the GNU General Public License. The 
Broyden mixing routine in PWscf is a simple module of 
about 20 statements that uses the blas and lapack li- 
braries. 



III. SPHERICAL DFT THEORY WITH 
MEAN-FIELD MIXING: TWO-COMPONENT 
UNITARY FERMI GAS 



y{m+l) _ y{m) 



(n) (^g^ As a first application of the standard Broyden's 

method in nuclear structure calculations, we describe the 
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numerical solution of a two-component superfluid in a 
spherically symmetric harmonic trap using the superfluid 
local density approximation (slda) of dft. The system 
is described by a three-parameter local energy density 
functional [l^ adjusted to results of Monte-Carlo cal- 
culations for homogeneous matter. The functional de- 
pends on the three local densities: p{r) (particle density), 
T(r) (kinetic energy density), and K(r) (pairing tensor), 
which are constructed from the usual Bogoliubov two- 
component quasi-particle wave functions [Uk{r), Vfe(r)]: 

PtM = E Pkl'MEk), Pi(r) ^ E \Vk\'M-E,), 

k k 

TT(r) = E |VC/feP/T(i?fe), nir) - ^ \^Vk? hi-Ek). 

k k 

-ir)^^U,V:^^<^^^^l^Jl^ (9) 

k 

where k runs over the quasi-particle states, p = Pt + 
Pl, and the function friE) = [I + exp{E/T)]-^ is the 
thermal distribution function for fermions. The resulting 
Kohn-Sham equations can be written in the usual hfb 
form, 

with the HFB matrix given by 

\ At(r) -h(r) + Hj 

where the particle-hole potential r(r), which enters the 
particle-hole Hamiltonian 

h{r) = -^V' + rir), (12) 

and the pairing potential A(r) are both functions of the 
local densities. (The value of the effective mass m* in 
Eq. has been determined in Ref. JLS]-) The intro- 
duction of two chemical potentials /i| and eschews 
the need to use a different formalism (Pauli blocking) for 
deahng with odd systems [13;, 14]. 

The spherical hfb equations pO|) are solved by dis- 
cretizing 7i in the radial coordinate r using the so-called 
DVR basis [3, [H, [13 • We use two sets of abscissa - 
{fo,n} {n = l,---,^r) for even-£ and {t"!,™} for odd-^ 
partial waves. Thus, to represent the functions A(r), 
r(r), etc., we need only store 2Nr sets of data. 

The chemical potentials for the two species (p| and 
Hi) can be viewed as the Lagrange multipliers that enter 
the grand thermodynamic potential £ — ~ fi^pi. In 
the asymmetric situation, Sp, = p.^ — pi ^ 0, the spectra 
for the species are shifted so that the thermal distribu- 
tion functions friE) appearing in ^ results in different 
numbers of each species. At zero temperature this pro- 
duces a quantized step as 6p approaches the pairing gap. 
By introducing a small temperature, however, this dis- 
continuity can be smoothed without noticeably affecting 



the physics. This smoothing is important in order for 
Broyden's method to achieve convergence. 

Instead of treating the chemical potentials as external 
parameters and adjusting them separately during the it- 
eration process, we have found that they may be adjusted 
during the iteration in such a way that the process con- 
verges to a solution with the desired particle number con- 
straints. After some experimenting, we have found that 
- in order to speed up the convergence - the chemical po- 
tentials should be updated after every iteration according 
to a simple ansatz: 

N'^ — N 

p.^p. + V V TF(A^g), (13) 

where a =| or f, is the desired particle number, and 
PTpiN) is some typical positive scale for the chemical 
potentials, and is some typical positive scale for the 
particle numbers (such as the total desired particle num- 
ber). The parameter 77 can be adjusted to ensure con- 
vergence. Note that the computed particle number No- 
enters in exactly one place so that the update is zero if 
and only if — iV°. Thermodynamics ensures that the 
chemical potential is adjusted in the appropriate direc- 
tion, so as long as the scales in Eq. are appropriately 
chosen the process converges nicely. The only problems 
we have encountered are: (i) If the updated steps are too 
large, the system might oscillate wildly. This is easily 
compensated for by choosing a smaller value of rj] (ii) If 
the temperature is too low, the resulting function may be 
discontinuous, in which case the update (|13p may cause 
an oscillation. Introducing a small, but finite, tempera- 
ture (smaller than any relevant scale) cures this problem. 

The Broyden vector V'-™-' in our slda calculations con- 
tains the values of r(r) and A(r) on two DVR grids, and 
two chemical potentials: 

V = {r(ro,„), r(ri^„), A(ro,„), A(ri^„), pi^p^}. (14) 

Thus, for a basis of size Nr, the size of the Broyden vector 
is = iNr + 2. 

We found that in order to ensure consistency in the 
above scheme, all parameters used in the iteration should 
appear in the Broyden process. Thus, even though our 
description in terms of the DVR basis is somewhat redun- 
dant with the function evaluated at two sets of abscissa 
{''o.n} 3'iid {ri „}, one must keep track of all the infor- 
mation describing the Jacobian. Also, one must be ex- 
tremely careful to include any external parameters that 
are adjusted during the iteration in the Jacobian, as well 
as any internal parameters that are allowed to relax, but 
which might add hysteresis to the system. The output of 
a single iteration must be a smooth and uniquely defined 
function of only the input vector Vj^^^ . 

Figure [1] displays some typical results. Here we solve 
the problem of 60, 59, and 45 particles in a harmonic trap. 
We consider the cases with equal numbers (N^ , N^) — 
(30, 30), one particle extra (A^t^ ^i) = (30, 29), and half- 
paired (A^T'^i) = (30,15). The size of the Broyden 
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FIG. 1: (Color online) Comparison between the linear mix- 
ing algorithm and the standard Broyden's method in SLDA 
for systems of (iVr.^-l) = (30,30) (top), {Nt;,Ni) = (30,29) 
(middle), and {N-f,Ni) = (30,15) (bottom). The basis size 
is Nr=^0 and the convergence parameter of Eq. (|13p is ri=l. 
For the asymmetric systems, the equations exhibit a disconti- 
nuity responsible for small plateaus in the plots. Introducing 
a small temperature T = 0.02eF restores the continuity. Once 
the approximate solution is found, the temperature may be 
set back to zero. 



vector is N—162. We see that in all cases, Broyden's 
method performs substantially better than the linear 
mixing method. 

Note that there are many cases in which the full-step 
algorithm (a=l) does not converge. This is typically be- 
cause the system exhibits some form of oscillation and 
the full step overshoots each time. In our problem this is 
due to the chemical potential update, and indicates that 
we should probably choose a smaller factor rj in Eq. p3p . 
Broyden's method, however, automatically detects this 
and compensates accordingly, thus achieving good con- 



vergence even if rj is not properly chosen. 

In the odd-particle case (30, 29), a plateau appears 
where the convergence improves very slowly. This marks 
the point at which the chemical potentials adjust so as 
to split one state off from the Fermi sea. Prior to this, 
the large-scale structure of the initial state is corrected, 
but once this occurs, the chemical potentials must be 
finely tuned so that exactly one extra state is occupied. 
Broyden's method has the same problem because the ap- 
proximate gradient changes rapidly. To deal with this, a 
small temperature is assumed. This makes the functions 
smooth and facilitates convergence. Once the steps are 
small enough that the energies do not cross zero, the os- 
cillations terminate and convergence is resumed. If prob- 
lems persist, the strategy is to start with a fairly large 
temperature to isolate the approximate values of the 
chemical potentials (only achieving rough convergence), 
and then iterating with these values at much lower tem- 
peratures to reach the desired tolerance. 



IV. SYMMETRY-UNCONSTRAINED HFB 
THEORY WITH MEAN-FIELD MIXING 

In this section we discuss the adaptation and per- 
formance of the modified Broyden's method in self- 
consistent nuclear structure calculations with the hfodd 
code of Ref. [H, [ll| . The code solves the Skyrme-HFB 
problem by expanding the solutions on the anisotropic 
cartesian harmonic oscillator (ho) basis, hfodd does 
not have any self-consistent symmetry built in; hence, 
all deformation degrees of freedom can be taken into ac- 
count. As the time- reversal symmetry can be broken in 
HFODD, the code can be used to describe nuclear rotation 
using the cranking approximation. 

The building blocks of the Skyrme-HFB approach are 
the particle density p(r, r') and the pairing density 
p{r,r'). The latter is related to the pairing tensor k 
of Sec, mil bv an anti- unitary transformation (see, e.g., 
0). ^ 

From p and p, other local densities and currents are 
constructed by considering spin and isospin degrees of 
freedom, and by taking derivatives up to the second or- 
der @. By coupling the local densities and currents, one 
constructs the energy density £{r) that is scalar, time- 
even, and isoscalar. In its current version, hfodd does 
not consider explicit proton-neutron mixing; hence, the 
HFB equations are solved for protons and neutrons sepa- 
rately. A variation of the total energy with respect to the 
single-particle wave functions yields the mean field h and 
pairing field h that enter the hfb Hamiltonian matrix: 



n = 



h{r) - A h{r) 
h{r) -h{r) + X 



(15) 



The particle-hole mean-field Hamiltonian can be ex- 
pressed as: 

h{r) = _ A+r^^+rg'^^'-fFf ™+rf'^-hC/'^™', (i6) 



2m 
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where the subscript t=0 (1) denotes isoscalar (isovector) 
fields, while the superscript even (odd) labels time-even 
(time-odd) fields. 

The pairing field h in hfodd is related to the usual 
pairing field A in Eq. by the same anti-unitary trans- 
formation that relates p and k. In our hfb calculations, 
we employed the density-dependent delta interaction in 
the particle-particle channel. The corresponding particle- 
particle mean-field Hamiltonian reads: 



Mr) = 



1 - Vi 



p{r) 



Po 



Pir), 



(17) 



where po^O.lGfm^'^. It this paper we adopted the value 
of Vi—0.5 corresponding to the so-called mixed pair- 
ing (see Ref. 21] and references therein). Due to the 
zero-range of the pairing interaction, a cut-off energy of 
60 MeV is used when summing up the contributions of the 
HFB quasi-particle states to the density matrices. For a 
given cut-off, the pairing strength Vb is determined such 
as to reproduce the experimental neutron pai ring gap 
in i20Sn, A„=1.245MeV. As discussed in Ref. [ij, this 
renormalization procedure gives very similar results as 
the pairing regularization method [23| . 

The mean-field Hamiltonian (fTB|) depends on the fol- 
lowing self-consistent potentials built of the local densi- 
ties and currents: 



-peve7i 
^ t 



(r) 



-V ■[Mt{r)V] + Ut{r) 



(18) 



+ - Va- Btir)+ Bt{r)- V ct 
Tf\r) =-V-[((T-C(r)V)]+cr.Si(r) 
+ l(Va./,(r) + 7,(r).V). 



The potentials Ut{r) and Mtr) are real scalar fields (1 
component), St('r), Ct{r) and It{r) are real vector fields 

(3 components), and B t{r) is a real second-rank tensor 
field (9 components). The pairing potential ht{r) is a 
two-component complex scalar field. 

The HO eigenfunctions are Hermite polynomials 
weighted by a Gaussian factor. Consequently, all inte- 
grations in HFODD are carried out by using the Gauss- 
Hermite quadrature pjj] . The mesh of the Gauss-Hermite 
nodes defines the discretized grid, and it is sufficient to 
define mean-field potentials on this grid to compute the 
HFB matrix. 

The Broyden vector V^f™) in hfodd contains the val- 
ues of the 44 self-consistent field components on the 
Gauss-Hermite grid rj: 

V = {Utin), Mt{ri),-Etiri), Ctin), Itin), 

lBtiri),^ht{ri)],^[ht{ri)]Y (19) 

The advantage of dealing with the fields is that one does 
not have to worry about the unitarity condition for the 
generalized density matrix. 



If Nx, Ny, and are the numbers of Gauss-Hermite 
nodes in the x-, y- and z-direction respectively, the total 
size of the Broyden vector in HFODD is = (N^NyNz) x 
44. For a spherical HO basis with Nose— ^"2, the number 
of Gauss-Hermite points required to perform exact inte- 
grations is Nx ~ Ny = Nz = 26, and the size of the 
Broyden vector becomes A^=773,344. In double preci- 
sion fioating point arithmetics, this represents approxi- 
mately 6.2 MB of memory. For M=8 vectors retained in 
the Broyden history, this represents a 100 MB memory 
overhead. Calculations involving larger bases are often 
required for studies of, e.g., fission pathways, and it is 
not unusual to deal with iVosc=20 for which the size of 
the Broyden vector becomes as large as iV=3,259,872, 
corresponding to ~416MB of additional memory. 

Figure [5] illustrates the use of Broyden's method for 
a cranked Hartree-Fock calculation (without pairing) for 
the lowest high-spin superdeformed band of ^^^Tb. The 
SKM* Skyrme functional 24] has been used. The rota- 
tion is generated by adding to the mean-field Hamilto- 
nian (jl6|) the cranking term —ujyjy, with ujy being the 
rotational frequency. As the time-reversal symmetry is 
broken in hf, both time-even and time-odd fields appear. 
The calculation at /ki;j,=0.5 MeV was initialized with the 
unconstrained hf results for the superdeformed vacuum 
configuration of ^^^Tb (with all the lowest single-particle 
HF routhians occupied) at LOy—0. The calculations at 
?iWj^=1.0MeV were initialized with the converged results 
at hiOy—O.^yieW . The same stretched HO basis contain- 
ing A'osc=15 shells and a deformation of /32=0.61 and 
Pi— 0.1 was used in both cases. The advantage of the 
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FIG. 2: (Color online) Comparison between linear mixing 
(q:=0.5, black symbols) and the modified Broyden's method 
(a=0.7, M—7) in hfodd for a high-spin superdeformed band 
in ^^^Tb at two values of rotational frequencies huoy^Q.b MeV 
and 1 MeV. 

modified Broyden's method over the simple mixing is ev- 
ident. While in the standard calculations the required 
tolerance of 10~^ is achieved after TO'^30 iterations, us- 
ing the linear mixing, the same outcome is achieved after 
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130 steps. 

Figure[3]illustrates the constrained hfodd calculations 
with pairing for the superdeformed minimum in ^^^Dy 
at the quadrupole moment of Q2=20eb (using the SLy4 
Skyrme functional [1^) and also for the spherical ground- 
state of ^°*Pb (using the SkP Skyrme functional [2^). 
Both calculations were carried out with mixed pairing 
and the Lipkin-Noeami (ln) particle-number corrections 
(see, e.g., Ref. [13, H^). In order to obtain stability of the 
results, the proton and neutron A2 parameters of the ln 
model had to be incorporated into the Broyden vector. 
Again, as in Fig. [2l the modified Broyden's method pro- 
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FIG. 3; (Color online) Comparison between linear mix- 
ing {a— 0.5, black symbols) and Broyden's method {a=0.7, 
M—7) in HFODD for the superdeformed minimum in ^^^Dy 
at the quadrupole moment of Q2=20eb (SLy4 Skyrme func- 
tional, Nosc—U) and the spherical ground state of ^°*Pb (SkP 
Skyrme functional, No3c~12). 

vides an impressive improvement over the linear mixing 
scheme. 

Before concluding this section, we would like to make 
an important technical remark. In cartesian coordinates, 
the HO wave-functions are products of the Hermite poly- 
nomials and the Gaussians. In the code hfodd the Gaus- 
sians do not have to be calculated explicitly, as they are 
already incorporated into the Gauss-Hermite quadrature 
weights [31 . However, this implies that the mean fields 
used in hfodd are Gaussian-scaled. Consequently, they 
can take non-zero (and possibly large) values at some 
points of the Gauss-Hermite grid where the physical po- 
tentials practically vanish. Using the scaled potentials 
as inputs to Broyden's method can be a source of very 
significant numerical instability. Indeed if the scaled po- 
tentials (or densities) are used, numerical and physically 
insignificant fluctuations of the real potentials are arti- 
ficially enhanced, thereby introducing considerable (and 
artificial) numerical noise. This instability is not present 
when the physical mean-field potentials are used. 



V. SYMMETRY-CONSTRAINED HFB THEORY 
WITH HAMILTONIAN MATRIX MIXING 

In this section we discuss the performance of the mod- 
ified Broyden's method as applied to the axial hfb prob- 
lem. To this end, we employ the hfbtho code [H, |2§| 
that solves the Skyrme-HFB equations by expanding the 
eigenstates in the cylindrical transformed harmonic oscil- 
lator (tho) basis. In the examples shown in this section, 
we use a SLy4 Skyrme energy density functional aug- 
mented by mixed pairing p7p . The hfb Hamiltonian 
matrix ([T5)) is represented in a tho basis. 

The Broyden vector \/(™) in hfbtho contains the ma- 
trix elements of the self-consistent fields h and h for neu- 
trons and protons: 

V^{h-,h^^,h^^,h^^}. (20) 

Since h and h are hermitian, i < j. The chemical po- 
tentials for neutrons and protons are adjusted at each 
step m to the correct number of particles. The matrix 
formulation (PO)) makes it easy to efficiently incorporate 
non-local terms that appear, e.g., in the particle-number 
restoration schemes. Some results of our Broyden hf- 
btho calculations have been reported in Ref. [30] . 

Since hfbtho strictly preserves axial symmetry, par- 
ity, and time reversal, the number of independent matrix 
elements is significantly reduced as compared to the full 
symmetry-unconstrained case of Sec. IIVI For instance, 
for the basis space of A'osc=20 principal oscillator shells, 
the length of the Broyden vector is iV=261,228. For 
the modified Broyden's method with M=8, the required 
memory increase in hfbtho as compared to the simple 
linear mixing algorithm is usually less than 33MB, i.e., 
is practically negligible. 

In the implementation of the linear mixing proce- 
dure in HFBTHO, the mixing parameter a varies. Ini- 
tially its value is set at a=0.1. In the next iteration, 
if the largest element of F^"^\ maxji^*^™^!, is less than 
maxji^*^™"^^ I, a is multiplied by a factor of 1.13 imtil it 
reaches the maximum value of q;=1.0. When the condi- 
tion max|i^(™^ |>max|i^'^™~^^ I is met, a is returned to its 
initial value of 0.1, and the process starts all over again. 
Such a technique has yielded stable results for all cases 
considered. 

Figure |3] compares linear mixing and Broyden's 
method for the spherical nucleus ^^"Sn, which is often 
used to determine the pairing strength of the functional. 
Both calculations are initiated from the same Woods- 
Saxon fields Usually convergence of all physical 
quantities of interest is achieved when max|i^(™^| is less 
than 10~^. For the nucleus ""^^"Sn, which represents a 
typical case, Broyden's method increases the convergence 
rate by a factor of '^3, as compared to the linear mixing. 
While the absolute convergence rate does depend on an 
individual nuclear configuration, Broyden's method has 
proved faster than the linear mixing in all the cases in- 
vestigated. In particular, for some pathological configu- 
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rations the improvement is staggering. Such an example 
is shown in Fig. 2] (bottom) for the prolate configuration 
in ^^"^Rn. Here, the efficiency of Broyden's algorithm is 
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FIG. 4: (Color online) Comparison between linear mixing 
(dotted line) and Broyden's method in hfbtho for ^^"Sn (top) 
and ^^''Rn (bottom). The largest element of F^™' is shown 
as a function of m. The modified Broyden mixing was carried 
out with a=0.7, 100=0.01, and Wn~l for M=3 (dashed line) 
and M=7 (solid line). The inset shows the slow convergence 
of the linear mixing for the pathological case of ^^■'Rn. Here, 
the tolerance fO~^ is reached after 4,345 iterations. 

almost a hundred times better than the linear mixing. 
An additional 30% gain can be made by increasing the 
number of vectors retained in the Broyden history from 
A/=3 to M==7. 

Broyden's method proves especially helpful for con- 
strained calculations in which the total energy is cal- 
culated as a function of some collective variables (usu- 
ally expectation values of selected one-body operators). 
A typical example is shown in Fig. [5] that displays 
the quadrupole deformation energy curve for ^^^Ra. 
The quadrupole deformation (3 is related to the total 
quadrupole moment Q20 via the usual relation: 



Q2 



(21) 



Starting from the converged solution at the neighbor- 
ing point of the deformation curve, the calculation con- 
verges in no more than 20-30 iterations. The number 



of iterations is reduced by a factor of 2-3 with respect 
to the linear mixing used under the same conditions. 
Since constrained calculations involve a large number of 
points along the collective path, and in many applica- 
tions several collective coordinates are involved, such a 
gain greatly improves the scalability of computations. 
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FIG. 5: (Color online) Energy curve of ^^^Ra versus 
quadrupole deformation obtained in the constrained hfbtho 
calculations. The local minima (maxima) are marked by dots 
(squares). While in the unconstrained calculations the linear 
mixing converges at the local minima only, Broyden's method 
yields both local minima and maxima, as indicated by arrows. 

A feature of Broyden's method, which is not present 
when using the linear mixing, is the ability to converge 
at both minima and maxima of the energy curve, which 
represent solutions of unconstrained hfb calculations. 
The extrema in the energy curve of ^^^Ra are marked 
in Fig. [5] by black symbols. Starting the unconstrained 
calculation from a constrained solution at (3=0.25 and us- 
ing the linear mixing, one obtains the local minimum at 
/3«0.05. Under the same conditions, the Broyden mix- 
ing yields the maximum at A similar situation 
happens when starting with an initial guess at zero de- 
formation: the linear mixing converges to the slightly 
deformed prolate minimum, while the Broyden mixing 
yields the spherical maximum of the energy curve. 

This interesting feature of Broyden's method may be 
particularly useful for many-body tunneling calculations 
(e.g., fission), since it allows to compute the potential 
barrier height with a much better accuracy than when us- 
ing the linear mixing in a constrained mode. However, it 
also implies that a more detailed knowledge of the poten- 
tial energy surface is required to compute unconstrained 
minima, since some of the converged solutions may turn 
out to be local maxima. A possible solution to this prob- 
lem is to make use of the curvature condition 



(m-1) 



y(m)\ g(m)_p(m)^Q (22) 



which guarantees that the Broyden step is downhill, i.e.. 
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toward the minimum. In practical calculations, we ac- 
cept the Broyden step if Eq. ([^ is satisfied. Otherwise, 
we apply the linear mixing while keeping the current iter- 
ation in the Broyden history. Such a procedure bypasses 
maxima or inflection points and properly converges to 
the nearest local minimum. 



VI. COUPLED-CLUSTER APPLICATIONS 

In this section we illustrate the efficiency of the mod- 
ified Broyden's method in the context of the nuclear 
CC theory. Coupled-Cluster theory originated in nu- 
clear structure; it was pioneered by Coester and Kiimmel 
[Sll,!!^! in the early 1960s. In the last two decades, it has 
been applied mostly in quantum chemistry and today 
it defines the state-of-the-art many-body theory in this 
field (see fs^l for a recent overview). Recently, the CC 
approach has seen a revival in nuclear structure [33 |. In 
[33] it was applied to the description of loosely bound and 
unbound helium isotopes; in [36l| it was shown that cou- 
pled cluster meets few-body benchmarks, and converged 
results for the ground state energy of '^"Ca were given; in 
coupled-cluster theory was extended to three-body 
Hamiltonians, and the first CC calculations with three- 
nucleon forces were presented. 

Within the CC theory, the nuclear ground-state wave 
function is written as 

|V)=e*|0), (23) 

where \4>) = OiLi ^1|0) is a single-particle product state 
expressed in some convenient representation and 

f = fi + fa + . . . + (24) 

is a correlation operator that can be expressed in terms 
of fc-body operators 

ii,...,Zfc;ai,...,afc 

(25) 

representing particle-hole (p-h) excitations with respect 
to the reference state |(/)). In Eq. ((24)) and in the follow- 
ing, fc, . . . label hole orbitals, while a,b,c, . . . refer to 
particle states. 

The many-body correlations atop the reference state 
are introduced through the exponentiated correlation op- 
erator. If no truncations are introduced in T, the theory 
is exact and the choice of starting reference state is arbi- 
trary. However, in practice, one typically restricts expan- 
sion ([23]) to the two leading terms, f w fi -f fa. Within 
the resulting coupled-cluster singles and doubles (ccsd) 
model, the CC equations can be written as 

E = (<^|i^0), (26) 
= mH\j>), (27) 
= (C'liJl^). (28) 



Here lipi^ 'T") = • ■ ■ ^Ii'^ii ■ ■ ■ o-i^ is a — nh ex- 
citation of the reference state {(p), and 

H = e"*i7e* = {^e^)^ (29) 

is the similarity-transformed Hamiltonian (note that H 
is non-Hermitian) . The last expression on the right-hand 
side of Eq. ([29l) indicates that only fully connected dia- 
grams contribute to the construction. This ensures that 
no unlinked diagrams enter the CC wave function, regard- 
less of the truncation level in T, and this property makes 
CC theory size extensive [3S| ■ 

In order to determine the ground state energy E in 
Eq. ([^ . the p-h excitation amplitudes t"; and tfj have 
to be evaluated through a nonlinear set of equations ([TT)) 
and (|28p . By making use of diagram rules and defin- 
ing intermediates by writing the diagrams in factorized 
form, the equations for the amplitudes can be written in 
a quasi-linearized form. As an example, Eq. (|27p can be 
written as: 

n _i r '^f* h'- +™ o.ie j-m _, Ze ^mi _ }_Zie ^mn 

Ja'^a^e ""m^a "ma^e ' "'m''ea ^ mn'^ae 

+ l<LtlJ. (30) 

In ([50]) . /* is the Fock matrix and I^, hl^, and are 
the intermediates [38| . The non-linearity is hidden in the 
intermediates which depend on the amplitudes tf and f^j. 

Typically the number of unknowns is too large to allow 
for a direct solution of the CC equations; hence, one has to 
resort to iterative approaches. In the largest calculation 
we have performed to date fsT], we solved Eqs. (|27|) and 
(PS)) for ~ 10® number of unknowns, on 1000 2GB parallel 
processors. Each iterative step costs, at the CCSD level, 
~ n^rip computational cycles [nh/rip is the number of 
holc/particle orbitals). 

Considering the number of unknowns and the highly 
non-linear nature of CC equations, the use of iterative 
convergence accelerators, such as the modified Broyden's 
method, is a must. In our nuclear CC calculations, the 
Broyden vector is given by the amplitudes: 

v^{t:i,t^]. (31) 

Figure [5] compares the performance of the modified 
Broyden's method with M=4, 8, and 12 to the simple 
mixing ansatz To solve the CCSD equations for the 
CCSD ground state energies of ^®Si (top) and '^^Ni (bot- 
tom), we have used an oscillator representation of the 
Hamiltonian in a model space of iVosc=10 major oscilla- 
tor shells with /ia;osc=32 MeV. The Hamiltonian consists 
of the kinetic-energy operator (minus the center-of-mass 
correction) and the N^LO nuclcon-nuclcon interaction of 
Ref. [3^. It is seen that utilizing Broyden's method in- 
creases the convergence dramatically as compared to the 
simple mixing. It is found that for M=12, iterations con- 
verge more rapidly, and the tolerance of 10~* is reached 
typically after 30-40 iterations. For the linear mixing to 
achieve similar accuracy, 130-140 iterations are required. 
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FIG. 6: (Color online) Convergence of CCSD ground state 
energy of ^®Si (top) and ^^Ni (bottom) using the modified 
Broyden's method with Af=4, 8, and 12 compared to the 
simple mixing with a=0.7. The calculations were carried out 
in a model space of A'osc^lO major oscillator shells. 



VII. CONCLUSIONS 

Broyden's method, a quasi-Newton method for the nu- 
merical solution of nonlinear equations in many vari- 
ables, is widely used in quantum chemistry to perform 
electronic-structure calculations. In this work, we ap- 
plied Broyden's method (both in its standard and modi- 
fied variants) to several nuclear structure problems. The 
examples range from several dft applications in various 
geometries to ab-initio CC calculations. 

The number of unknown variables (i.e., the size of the 
Broyden vector) in our calculations ranges from ^^200 in 
SLDA calculations with the standard Broyden's method, 
to ~300,000 in hfbtho, -3,000,000 in hfodd, and ~10^ 
in CC applications of the modified Broyden's algorithm. 



Much faster convergence has been achieved in compari- 
son with the linear-mixing procedure which is often used 
in such types of calculations. We note that Broyden's 
method often achieves superlinear convergence, although 
mildly. If the iteration is not properly structured, then it 
might not correspond to a contractive mapping and may 
not converge. By constructing an approximation to the 
Jacobian, Broyden's method can determine the downhill 
direction even if the iteration itself does not provide this 
information. 

The downhill direction chosen by Broyden's method 
is towards a fixed point, not necessarily towards a mini- 
mum of the functional. Thus, the iterative process might 
easily converge to a maximum or a saddle point even 
when the iteration would naturally head towards a min- 
imum. If the goal is to find the minimum, a combina- 
tion of the linear mixing and Broyden's method can be 
efficiently used, as demonstrated in the example of con- 
strained HFB calculations of Fig. [5] However, the above 
feature of Broyden's method to converge to maxima or 
saddle points can be taken advantage of if the theoretical 
objective is to determine local maxima or saddle points 
(e.g., when studying fission pathways). 

In summary, Broyden's method is easy to incorporate 
into fixed-point codes, and provides impressive perfor- 
mance improvements in nuclear many-body calculations. 
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